/*
 * File:   InitParams.cpp
 * Author: georgy
 *
 * Created on February, 2012
 */

#include <math.h>
#include "MDAlgorithms.h"



// single step of simulation algothim.
void simulation_step(std::vector<fltParticle>  &p,
                     float                     dt,     // time step.
                     float                     cutOff, // cut-off.
                     float                     L)      // Linear size of cube.
{
    static float scale = 0.9;

    for (int i = 0; i < p.size(); ++i) { // Scale down velocities.
        p[i].vx *= scale;
        p[i].vy *= scale;
        p[i].vz *= scale;
    }

    // Verlet algorithm.
    for (int i = 0; i < p.size(); ++i) {
        p[i].x += (p[i].vx * dt) + (p[i].ax * dt * dt * 0.5);
        p[i].y += (p[i].vy * dt) + (p[i].ay * dt * dt * 0.5);
        p[i].z += (p[i].vz * dt) + (p[i].az * dt * dt * 0.5);

        p[i].vx += (p[i].ax * dt * dt * 0.5);
        p[i].vy += (p[i].ay * dt * dt * 0.5);
        p[i].vz += (p[i].az * dt * dt * 0.5);
    }

    calc_accelerations(p, cutOff);

    for (int i = 0; i < p.size(); ++i) {
        p[i].vx += (p[i].ax * dt/2);
        p[i].vy += (p[i].ay * dt/2);
        p[i].vz += (p[i].az * dt/2);
    }

    // Check boundary conditions (if particle reaches the wall).
    for (int i = 0; i < p.size(); ++i) {
        if (p[i].x >= L) {
            p[i].x = L;
            p[i].vx = -p[i].vx;
        }

        if (p[i].x < 0.0) {
            p[i].x = 0.0;
            p[i].vx = -p[i].vx;
        }

        if (p[i].y >= L) {
            p[i].y = L;
            p[i].vy = -p[i].vy;
        }

        if (p[i].y < 0.0) {
            p[i].y = 0.0;
            p[i].vy = -p[i].vy;
        }

        if (p[i].z >= L) {
            p[i].z = L;
            p[i].vz = -p[i].vz;
        }

        if (p[i].z < 0.0) {
            p[i].z = 0.0;
            p[i].vz = -p[i].vz;
        }
    }
}



void calc_accelerations(std::vector<fltParticle>  &p,
                        float                     cutOff) // cut-off distance.
{
    float dx     = 0.0;
    float dy     = 0.0;
    float dz     = 0.0;
    float fx     = 0.0;
    float fy     = 0.0;
    float fz     = 0.0;
    float r2     = 0.0;
    float attr   = 0.0;
    float repel  = 0.0;
    float fOverR = 0.0;

    // Colculate accelerations using Lennard-Jones potential.
    for (int i = 0; i < p.size(); ++i) {
        for (int j = 0; j < i; ++j) {
            dx = p[i].x - p[j].x;
            dy = p[i].y - p[j].y;
            dz = p[i].z - p[j].z;

            r2 = dx*dx + dy*dy + dz*dz;

            if (r2 < (cutOff * cutOff)) {
                attr = (1.0 / r2) * (1.0 / r2) * (1.0 / r2);
                repel = attr * attr;
                fOverR = 24.0 * ((1.3 * repel) - attr) * (1.0 / r2);
                fx = fOverR * dx;
                fy = fOverR * dy;
                fz = fOverR * dz;
                p[i].ax += fx;  // (mass = 1).
                p[i].ay += fy;
                p[i].az += fz;
                p[j].ax -= fx;  // Newton's 3rd law.
                p[j].ay -= fy;
                p[j].az -= fz;
            }
        }
    }
}



void initialize_conditions(std::vector<fltParticle>  &p,
                           float                     L,   // Linear size of cube.
                           float                     T)  // Temperature.
{
    // Set the particles in the nodes of lattice.
    int     N          = p.size();
    double  unitVolume = floor((L * L * L) / N);
    double  unitLength = floor(pow(unitVolume, 1/3.0));
    double  x          = unitLength / 2;
    double  y          = unitLength / 2;
    double  z          = unitLength / 2;

    for (int i = 0; i < N; ++i) {
        p[i].x = x;
        p[i].y = y;
        p[i].z = z;
        x += unitLength;

        if (x > L) {            // when right border of lattice is reached.
            x = unitLength / 2;
            y += unitLength;
        }

        if (y > L) {            // when bottom of lattice is reached.
            y = unitLength / 2;
            x = unitLength / 2;
            z += unitLength;
        }
    }

    // Normalize velocities.
    // Scale the velocity to satisfy the partition theorem.
    int     nv = 3 * N;
    int     ng = nv - 6;
    double  ek = 0;
    double  vs = 0;

    for (int i = 0; i < N; ++i) {
        ek += pow(p[i].vx, 2) + pow(p[i].vy, 2) + pow(p[i].vz, 2);
        vs = sqrt(ek * nv / (ng * T));

        for (int j = 0; j < N; ++j)
        {
            p[j].vx /= vs;
            p[j].vy /= vs;
            p[j].vz /= vs;
        }
    }
}



float calc_temperature(std::vector<fltParticle>  &p)
{
    double N = p.size();
    double K = 0.0; // Kinetic energy.
    double v = 0.0; // <v>.

    for (int i = 0; i < N; ++i) {
        v = sqrt(pow(p[i].vx, 2) + pow(p[i].vy, 2) + pow(p[i].vz, 2));
        K += 0.5 * v * v;  // (mass = 1) ;
    }

    // 3/2*NkT = <K>
    return (2 * K / (3 * N));
}

